ESTADISTICA PARA EL ANALISIS POLITICO II

Regresión Lineal Multivariada (I)

El Pipeline del Analísta Político

De manera abstracta podemos ver que Los analistas existen para brindar explicaciones y recomendaciones sobre algún asunto de interés. El decisor elije el curso de acción sabiendo que no puede esperar del analista información completa. El gestor se encarga de implementar las decisiones; esa implementación traerá nueva información y el analista vuelve a su trabajo.

Para acabar este curso sintiéndote que vas en camino a ser analista político, te recomiendo seguir estos pasos:

  1. Identificar algún patrón (no aleatorio) en la complejidad del mundo, que afectan los objetivos de las organizaciones: el problema. La identificación del problema depende de cómo entiendes al mundo, por lo que debes haber leído diversas teorías antes de aventurarte a identificar.

  2. Estructurar lo identificado (el problema) en preguntas de investigación.

  3. Revisar la literatura para saber qué respuestas hayan dado otros investigadores al problema planteado.

  4. Plantear explicaciones al problema: hipótesis.

  5. Organizar una estrategia para probar, o invalidar, las hipótesis.

  6. Aplicar la estrategia, entender y validar los resultados. Asegurar la replicabilidad de la estrategia.

  7. Elaborar síntesis de lo actuado; proponer explicaciones a lo encontrado; y enviar recomendaciones.

Hay muchas opciones para la estrategia señalada en el punto 5. La elección de las mismas dependerá de la preparación del analista. En este curso, te daremos conocimientos para que no desestimes el camino basado en datos y/o evidencias empíricas. De ahí que parte de tu proyecto para este curso está basado en desarrollar los puntos anteriores.

Para mayores detalles, revisa el SILABO presente en Paideia y el campus virtual PUCP.

En esta sesión deseo que entiendas por qué es necesario ir más allá de la correlación (Pearson, Spearman, etc.) o las diferencias de valores centrales (t test, kruska wallis, etc.). Esta necesidad de ir más allá trae una técnica conocida como la regresión.

La sesión va en dos partes:

I. De Correlación a regresión

Trabajemos con estos datos:

library(rio)
topCloud='https://github.com/PoliticayGobiernoPUCP/'
courseCloud='estadistica_anapol2/raw/master/DATA/'
fileCloud='hsb_ok.xlsx'
hsb=import(paste0(topCloud,courseCloud,fileCloud))

Antes de correr cualquier análisis estadístico, debes revisar como el tipo de datos que tu archivo trae es reconocido por el R:

str(hsb)
## 'data.frame':    600 obs. of  15 variables:
##  $ ID    : chr  "1" "2" "3" "4" ...
##  $ SEX   : num  2 1 2 2 2 1 1 2 1 2 ...
##  $ RACE  : num  2 2 2 2 2 2 2 2 2 2 ...
##  $ SES   : num  1 1 1 2 2 2 1 1 2 1 ...
##  $ SCTYP : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ HSP   : num  3 2 2 3 3 2 1 1 1 1 ...
##  $ LOCUS : num  0.29 -0.42 0.71 0.06 0.22 0.46 0.44 0.68 0.06 0.05 ...
##  $ CONCPT: num  0.88 0.03 0.03 0.03 -0.28 0.03 -0.47 0.25 0.56 0.15 ...
##  $ MOT   : num  0.67 0.33 0.67 0 0 0 0.33 1 0.33 1 ...
##  $ CAR   : num  10 2 9 15 1 11 10 9 9 11 ...
##  $ RDG   : num  33.6 46.9 41.6 38.9 36.3 49.5 62.7 44.2 46.9 44.2 ...
##  $ WRTG  : num  43.7 35.9 59.3 41.1 48.9 46.3 64.5 51.5 41.1 49.5 ...
##  $ MATH  : num  40.2 41.9 41.9 32.7 39.5 46.2 48 36.9 45.3 40.5 ...
##  $ SCI   : num  39 36.3 44.4 41.7 41.7 41.7 63.4 49.8 47.1 39 ...
##  $ CIV   : num  40.6 45.6 45.6 40.6 45.6 35.6 55.6 55.6 55.6 50.6 ...

Casi todo salió numérico, pero no sabremos qué ajustar si no leemos el manual metodológico o el diccionario de datos o la metadata disponible.

De ahi que debemos pre procesar:

categoricals=c("SEX","RACE","SES","SCTYP","HSP","CAR")

hsb[,categoricals]=lapply(hsb[,categoricals], as.factor)

# nominales
hsb$SEX=factor(hsb$SEX,
                levels=c(1,2),
                labels=c("Male","Female"))

hsb$RACE=factor(hsb$RACE,
                levels=c(1,2,3,4),
                labels=c("Hispanic","Asian","Black","White"))

hsb$HSP=factor(hsb$HSP,
                levels=c(1,2,3),labels=c("General","Academic","Vocational"))

hsb$SCTYP=factor(hsb$SCTYP,
                levels=c(1,2),
                labels=c("Public","Private"))
# a ordinal:
hsb$SES=ordered(hsb$SES,
                levels=c(1,2,3),
                labels=c("Low","Medium","High" ))

La variable de interés:

Asumamos que nuestra variable de interés es el desempeño en matemáticas; así, nuestra variable dependiente está representada por la variable MATH.

A partir de ahí, consideremos que nos interesa saber la posible relación que pueda tener la variable que ha medido el desempeño en escritura; así, una variable independiente sería la representada por la variable WRTG. Hasta ahora sabemos que como son dos variables de tipo numérico debemos usar una correlación. La gráfica de correlación es esta:

library(ggplot2)

base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
scatter = base + geom_point()
scatter

Vemos que hay una aparente relación. Calculemos los indices de correlación:

f1=formula(~MATH + WRTG)

# camino parametrico
pearsonf1=cor.test(f1,data=hsb)[c('estimate','p.value')]

# camino no parametrico
spearmanf1=cor.test(f1,data=hsb,method='spearman')[c('estimate','p.value')]

Asumiendo un camino paramétrico, podemos pedir el coeficiente de Pearson, el cuál al ser calculado obtenemos 0.6326664 (con p-value= 0).

Si hubieramos seguido una ruta no paramétrica, informaríamos el coeficiente de Spearman:0.6415126 (con p-value= 0).

Ahora, consideremos que nos interesa saber a la vez la posible relación que pueda tener la variable que ha medido el desempeño en ciencias; así otra variable independiente sería la representada por la variable SCI. Como es otra variable numérica no podemos calcular la correlación de tres variables, pero podemos tratar de verlo visualmente:

  • En tres dimensiones:
library(scatterplot3d)
scatterplot3d(hsb[,c('SCI','WRTG','MATH')])

  • En dos dimensiones:
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SCI))

Calculemos las correlaciones:

f2=formula(~MATH+SCI)

# camino parametrico
pearsonf2=cor.test(f2,data=hsb)[c('estimate','p.value')]

# camino no parametrico
spearmanf2=cor.test(f2,data=hsb, method='spearman')[c('estimate','p.value')]

Podríamos calcular la correlación de SCI con MATH, obteniendo el Pearson (0.6495261, p-value= 0) y el Spearman (0.6551515,p-value= 0).

Visualmente vemos relación, pero no tenemos un coeficiente para medir ello.

Finalmente, ¿si quisiéramos ver si el sexo tiene algun rol en todo esto? Como ésta es una variable categórica y dicotómica, lo primero que puede venir a la mente es esta gráfica:

base=ggplot(data=hsb, aes(x=as.factor(SEX), y=MATH))
base + geom_boxplot(notch = T) +  geom_jitter(color="black", size=0.4, alpha=0.9)

Los boxplots tienen un notch flanqueando a la mediana, para sugerir igualdad de medianas si éstos se intersectan; de ahi que parece no haber diferencia sustantiva entre hombres y mujeres en cuanto a su desempeño en MATH.

Una alternativa al boxplot seria las barras de error:

library(ggpubr)
ggerrorplot(data=hsb, x = "SEX", y = "MATH")

En este último caso, si las lineas (denotado por las barras de error de la media) se intersectan, sugeriria que los valores medios (en este caso la media) podrian ser iguales.

Verificar si hay o no igualdad entre distribuciones depende si las variables se distribuyen o no de manera normal:

library(ggplot2)
ggplot(hsb,aes(x=MATH)) + 
  geom_histogram(aes(y = ..density..),bins = 20, fill='green') +
  stat_function(fun = dnorm, colour = "red",
                      args = list(mean = mean(hsb$MATH, na.rm = TRUE),
                                 sd = sd(hsb$MATH, na.rm = TRUE))) + 
  facet_grid(~SEX) + 
  coord_flip()

Nota que los histogramas de la data real tienen encima la curva normal que idealmente tendría esa data. La lejanía entre ellos, sugeriría no normalidad.

Se suele usar un qqplot para explorar la presencia/ausencia de normalidad:

# se sugiere normalidad si los puntos no se alejan de la diagonal.

library(ggpubr)
ggqqplot(data=hsb,x="MATH") + facet_grid(. ~ SEX)

Como ello no es fácil de discernir visualmente, tenemos por costumbre calcular algun coeficiente, como el Shapiro-Wilk:

library(knitr)
library(magrittr)
library(kableExtra)
f4=formula(MATH~SEX)


tablag= aggregate(f4, hsb,
          FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})

# para que se vea mejor:


shapiroTest=as.data.frame(tablag[,2])
names(shapiroTest)=c("W","Prob")


kable(cbind(tablag[1],shapiroTest))%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F, position = "left")
SEX W Prob
Male 0.9837903 0.0034565
Female 0.9790040 0.0001031

Esto nos sugiere un camino no paramétrico para ver la diferencia de valores medios, por lo que deberiamos usar la prueba de Mann-Whitney en vez de la prueba t para testaer la relación entre ambas.

tf4=t.test(f4,data=hsb)[c('estimate','p.value')]
wilcoxf4=wilcox.test(f4,data=hsb)['p.value']

La prueba no paramétrica no rechazaría la igualdad de valores medios (Mann-Whitney con p valor = 0.3085543).

  • Veamos como representar al sexo en nuestra gráfica entre WRTG y MATH:
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SEX))

  • Veamos como representar al sexo en nuestra gráfica entre WRTG, SCI y MATH:
base + geom_point(aes(size = SCI, color=SEX))

Otra alternativa puede ser:

base + geom_point(aes(color = SCI)) + facet_grid(~SEX)

Y claro:

paleta <- c("coral1","cyan" )
colors <- paleta[as.numeric(hsb$SEX)]
scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)

En todos los gráficos vemos que los hombres y las mujeres están distribuidos por todo el gráfico, lo cual nos sugiere que no hay diferencias aun en dimensiones mayores a dos. Sin embargo, no tenemos una medida de cuanto cada uno afecta a nuestra dependiente.

De ahi que necesitamos la regresión.

Regresión Lineal

La regresión es una técnica donde hay que definir una variable dependiente y una o más independientes. Las independientes pueden tener rol predictor, dependiendo del diseño de investigación, aunque por defecto tiene un rol explicativo.

La regresión sí quiere informar cuánto una variable (independiente) puede explicar la variación de otra (dependiente), de ahí que es una técnica para probar hipótesis direccionales o asimétricas (las correlaciones tiene hipótesis simétricas).

La regresión busca proponer un modelo, es decir una ecuación, que recoja como una variable explicaría a otra.Para nuestro caso la variable dependiente es MATH, y tendriamos hasta ahora tres modelos:

modelo1=formula(MATH~WRTG)
modelo2=formula(MATH ~ WRTG + SCI)
modelo3= formula(MATH ~ WRTG + SCI + SEX)

Por ejemplo, para la hipótesis ‘el nivel de desempeño en escritura afecta el desempeño en matemáticas’, la regresión arrojaría este resultado:


O en mejor version con ayuda de stargazer:

library(stargazer)
reg1=lm(modelo1,data=hsb)
stargazer(reg1,type = "html",intercept.bottom = FALSE)
Dependent variable:
MATH
Constant 19.769***
(1.633)
WRTG 0.612***
(0.031)
Observations 600
R2 0.400
Adjusted R2 0.399
Residual Std. Error 7.297 (df = 598)
F Statistic 399.110*** (df = 1; 598)
Note: p<0.1; p<0.05; p<0.01


Aquí ya sabemos algo interesante, primero que WRTG tiene efecto, pues es significativo (indicado por los dos asteriscos); segundo, que ese efecto es directo, pues el coeficiente calculado es positivo; y tercero que la magnitud de ese efecto es 0.612, lo que indica cuanto aumenta MATH en promedio cuando WRTG se incremente en una unidad.

Esto es información suficiente para representar esa relación con una ecuación. Como la ecuación sólo tiene una variable independiente podemos producir una recta sobre el gráfico de correlación:

ggplot(hsb, aes(x=WRTG, y=MATH)) + 
  geom_point()+
  geom_smooth(method=lm)

Esa recta podemos representarla así:

\[ MATH= 19.7690323 + 0.6123904 \cdot WRTG + \epsilon\]

El Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado, de ahi la presencia del \(\epsilon\). Justamente el R cuadrado ajustado (0.4002668) nos brinda un porcentaje (multiplicalo por 100) que da una pista de nuestra cercanía a una situación perfecta (cuando vale 1).

Y sí queremos ver el efecto de SCI?

reg2=lm(modelo2,data=hsb)
stargazer(reg2,type = "html",intercept.bottom = FALSE)
Dependent variable:
MATH
Constant 10.629***
(1.630)
WRTG 0.377***
(0.033)
SCI 0.415***
(0.033)
Observations 600
R2 0.524
Adjusted R2 0.523
Residual Std. Error 6.505 (df = 597)
F Statistic 328.846*** (df = 2; 597)
Note: p<0.1; p<0.05; p<0.01

En este caso, la regresión tendrá una formula con dos variables explicando la dependiente, así que en vez de producir una línea buscará producir un plano:

library(scatterplot3d)
G  <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)

Este plano podemos representarlo así:

\[ MATH= 10.6285904 + 0.3765304 \cdot WRTG + 0.4152733 \cdot SCI + \epsilon\]

Nuevamente, el Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado en forma de plano. De igual manera el R cuadrado ajustado (0.5241861) nos da una pista de nuestra lejanía a una situación perfecta.

Es clave darse cuenta de otro detalle, que el coeficiente de WRTG ha variado en la fórmula ahora que está presente SCI ¿Por qué sucede esto? Veamoslo así: en el primer caso, WRTG y \(\epsilon\) buscaban representar la variabilidad en MATH, y ahora, en el segundo caso, viene SCI para mejorar esa explicación; es así que el peso de la explicación ahora se recalcula y el coeficiente de WRTG deja de explicar lo que le corresponde a SCI, y \(\epsilon\) también le entrega algo a SCI.

Como \(\epsilon\) no tiene coeficiente, representamos su variación usando el error típico de los residuos o residual standard error (RSE). Nótese que éste ha variado de un modelo ha otro, ahora tenemos un RSE menor. Aquí vale la pena preguntarse si esta disminución del error es significativa, obteniendo:


tanova=anova(reg1,reg2)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
Table de Análisis de Varianza
Res.Df RSS Df Sum of Sq F Pr(> F)
1 598 31,842.070
2 597 25,262.730 1 6,579.335 155.481 0


La comparación de modelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasar de un modelo a otro). Como la comparación es significativa (vea el P), rechazamos igualdad de modelos: el modelo 2 sí reduce el error al incluir una variable más.

Finalmente, veamos el rol de sexo:

reg3=lm(modelo3,data=hsb)
stargazer(reg3,type = "html",intercept.bottom = FALSE)
Dependent variable:
MATH
Constant 11.306***
(1.630)
WRTG 0.424***
(0.036)
SCI 0.375***
(0.035)
SEXFemale -1.922***
(0.582)
Observations 600
R2 0.533
Adjusted R2 0.530
Residual Std. Error 6.452 (df = 596)
F Statistic 226.510*** (df = 3; 596)
Note: p<0.1; p<0.05; p<0.01

Aunque no podemos graficar cuatro coordenadas, podemos usar elementos visuales:

library(scatterplot3d)
colors <- paleta[as.numeric(hsb$SEX)]
G  <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)

Nuestra nueva ecuación sería:

\[ MATH= 11.306349 + 0.4235786 \cdot WRTG + 0.3748025 \cdot SCI + -1.9219692 \cdot SEX + \epsilon\]

Nuevamente podemos ver si añadir SEXO en este modelo representa una mejora al anterior:

tanova=anova(reg1,reg2,reg3)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
Table de Análisis de Varianza
Res.Df RSS Df Sum of Sq F Pr(> F)
1 598 31,842.070
2 597 25,262.730 1 6,579.335 158.063 0
3 596 24,808.420 1 454.308 10.914 0.001



Finalmente, podemos resumir todo en esta tabla:


library(stargazer)
stargazer(reg1,reg2,reg3, type = "html", title = "Modelos planteadas",digits = 2, single.row = F,no.space = F,intercept.bottom = FALSE,
          dep.var.caption="Variable dependiente:",
          dep.var.labels="Desempeño en Matemáticas",
          covariate.labels=c("Constante","Desempeño en Escritura","Desempeño en Ciencias","SEXO (mujer)"),
          keep.stat = c("n","adj.rsq","ser"),df = F,
          notes.label = "Notas:")
Modelos planteadas
Variable dependiente:
Desempeño en Matemáticas
(1) (2) (3)
Constante 19.77*** 10.63*** 11.31***
(1.63) (1.63) (1.63)
Desempeño en Escritura 0.61*** 0.38*** 0.42***
(0.03) (0.03) (0.04)
Desempeño en Ciencias 0.42*** 0.37***
(0.03) (0.04)
SEXO (mujer) -1.92***
(0.58)
Observations 600 600 600
Adjusted R2 0.40 0.52 0.53
Residual Std. Error 7.30 6.51 6.45
Notas: p<0.1; p<0.05; p<0.01

Gráficamente:

library(ggplot2)
library(sjPlot)
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
plot_models(reg1,reg2,reg3,vline.color = "grey",m.labels=c("Modelo 1","Modelo 2","Modelo 3"))

Como vemos, el propósito del último gráfico es mostrar que dado que ninguna de las variables (y sus intervalos de confianza) tocan el valor cero (que significa que la variable no tiene efecto en la dependiente).

al INICIO